#ifndef DW_MATHS_POLYNOMIAL_INL
#	define DW_MATHS_POLYNOMIAL_INL

namespace dw {
	namespace maths {

//==============================================================================
// struct Polynomial< 1 >
//==============================================================================

//------------------------------------------------------------------------------
// Construction & Destruction
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
Polynomial< 1 >::Polynomial()
:	a(0.0f)
,	b(0.0f)
{
	;
}

//------------------------------------------------------------------------------
// Commands
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
uint32 Polynomial< 1 >::solve(float32& x)
{
	dwAssert(a != 0.0f, "coefficients are invalid");

	x = -b / a;

	return 1;
}

//==============================================================================
// struct Polynomial< 2 >
//==============================================================================

//------------------------------------------------------------------------------
// Constants
//------------------------------------------------------------------------------

/*static*/ const float32 Polynomial< 2 >::DiscriminantEspilon = 1e-6f;

//------------------------------------------------------------------------------
// Construction & Destruction
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
Polynomial< 2 >::Polynomial()
:	a(0.0f)
,	b(0.0f)
,	c(0.0f)
{
	;
}

//------------------------------------------------------------------------------
uint32 Polynomial< 2 >::solve(float32& x0, float32& x1)
{
	if (std::abs(a) < std::numeric_limits< float32 >::epsilon())
	{
		Polynomial< 1 > p1;

		p1.a = b;
		p1.b = c;

		return p1.solve(x0);
	}

	float64 delta =  b * b - 4.0 * a * c;

	if (delta < 0.0)
	{
		return 0;
	}

	float64 invHalfA = 0.5 / a;

	if (std::abs(delta) < DiscriminantEspilon)
	{
		x0 = (float32)(-b * invHalfA);

		return 1;
	}

	float64 sqrtDelta = std::sqrt(delta);

	if (invHalfA > 0.0f)
	{
		x0 = (float32)((-b - sqrtDelta) * invHalfA);
		x1 = (float32)((-b + sqrtDelta) * invHalfA);
	}
	else
	{
		x0 = (float32)((-b + sqrtDelta) * invHalfA);
		x1 = (float32)((-b - sqrtDelta) * invHalfA);
	}

	return 2;
}

//==============================================================================
// struct Polynomial< 3 >
//==============================================================================

//------------------------------------------------------------------------------
// Construction & Destruction
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
Polynomial< 3 >::Polynomial()
:	a(0.0f)
,	b(0.0f)
,	c(0.0f)
,	d(0.0f)
{
	;
}

//------------------------------------------------------------------------------
uint32 Polynomial< 3 >::solve(float32& x0, float32& x1, float32& x2)
{
	// if a ~= 0 then we can solve a Polynomial< 2 >
	if (std::abs(a) < std::numeric_limits< float32 >::epsilon())
	{
		Polynomial< 2 > p2;

		p2.a = b;
		p2.b = c;
		p2.c = d;

		return p2.solve(x0, x1);
	}

	//! @todo : check we can do this !!
	// if d ~= 0 then we can solve a Polynomial< 2 >
	/*
	if (std::abs(d) < std::numeric_limits< float32 >::epsilon())
	{
		Polynomial< 2 > p2;

		p2.a = a;
		p2.b = b;
		p2.c = c;

		return p2.solve(x0, x1);
	}
	*/

	float64 Q, R, Q3, R2, sQ, D, an, theta;
	float64 A2, a1, a2, a3;

	if (a != 1.0f)
	{
		a1 = b / a;
		a2 = c / a;
		a3 = d / a;
	}
	else
	{
		a1 = b;
		a2 = c;
		a3 = d;
	}

	A2 = a1 * a1;

	Q = (A2 - 3.0f * a2) / 9.0;

	R = (a1 * (A2 - 4.5 * a2) + 13.5 * a3) / 27.0;

	Q3 = Q * Q * Q;

	R2 = R * R;

	D = Q3 - R2;

	an = a1 / 3.0;

	if (D >= 0.0)
	{
		/* Three real roots. */

		D = R / std::sqrt(Q3);

		theta = std::acos(D) / 3.0;

		sQ = -2.0 * std::sqrt(Q);

		float64 tmpX0 = sQ * std::cos(theta) - an;
		float64 tmpX1 = sQ * std::cos(theta + two_pi_3) - an;
		float64 tmpX2 = sQ * std::cos(theta + two_pi_43) - an;

		x0 = (float32)std::min(std::min(tmpX0, tmpX1), tmpX2);
		x2 = (float32)std::max(std::max(tmpX0, tmpX1), tmpX2);

		if (tmpX1 < tmpX0 && tmpX0 < tmpX2)
		{
			x1 = (float32)tmpX0;
		}
		else if (tmpX0 < tmpX1 && tmpX1 < tmpX2)
		{
			x1 = (float32)tmpX1;
		}
		else
		{
			x1 = (float32)tmpX2;
		}

		return 3;
	}
	else
	{
		sQ = std::pow(std::sqrt(R2 - Q3) + std::abs(R), 1.0 / 3.0);

		if (R < 0.0)
		{
			x0 = (float32)((sQ + Q / sQ) - an);
		}
		else
		{
			x0 = (float32)(-(sQ + Q / sQ) - an);
		}

		return 1;
	}
}

//==============================================================================
// struct Polynomial< 4 >
//==============================================================================

//------------------------------------------------------------------------------
// Construction & Destruction
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
Polynomial< 4 >::Polynomial()
:	a(0.0f)
,	b(0.0f)
,	c(0.0f)
,	d(0.0f)
,	e(0.0f)
{
	;
}

//------------------------------------------------------------------------------
uint32 Polynomial< 4 >::solve(float32& x0, float32& x1, float32& x2, float32& x3)
{
	// if a ~= 0 then we can solve a Polynomial< 3 >
	if (std::abs(a) < std::numeric_limits< float32 >::epsilon())
	{
		Polynomial< 3 > p3;

		p3.a = b;
		p3.b = c;
		p3.c = d;
		p3.d = e;

		return p3.solve(x0, x1, x2);
	}

	//! @todo : check we can do this !!
	if (std::abs(d) < std::numeric_limits< float32 >::epsilon())
	{
		Polynomial< 3 > p3;

		p3.a = a;
		p3.b = b;
		p3.c = c;
		p3.d = d;

		return p3.solve(x0, x1, x2);
	}

	float32 roots[3];
	float64 c12, z, p, q, q1, q2, r, d1, d2;
	float64 c0, c1, c2, c3, c4;
	int i;

	c0 = a;

	// See if the constant term has vanished 
	// Make sure the quartic has a leading coefficient of 1.0 
	if (c0 != 1.0f)
	{
		c1 = b / c0;
		c2 = c / c0;
		c3 = d / c0;
		c4 = e / c0;
	}
	else
	{
		c1 = b;
		c2 = c;
		c3 = d;
		c4 = e;
	}

	// Compute the cubic resolvant 
	c12 = c1 * c1;
	p = -0.375 * c12 + c2;
	q = 0.125 * c12 * c1 - 0.5 * c1 * c2 + c3;
	r = -0.01171875 * c12 * c12 + 0.0625 * c12 * c2 - 0.25 * c1 * c3 + c4;

	Polynomial< 3 > cubic;

	cubic.a = 1.0f;
	cubic.b = (float32)(-0.5 * p);
	cubic.c = (float32)(-r);
	cubic.d = (float32)(0.5 * r * p - 0.125 * q * q);

	i = cubic.solve(roots[0], roots[1], roots[2]);

	if (i > 0)
		z = roots[0];
	else
		return 0;

	d1 = 2.0 * z - p;

	if (d1 < 0.0)
	{
		if (d1 > -std::numeric_limits< float64 >::epsilon())
			d1 = 0.0;
		else
			return 0;
	}

	if (d1 < std::numeric_limits< float64 >::epsilon())
	{
		d2 = z * z - r;
		if (d2 < 0.0)
			return 0;
		d2 = std::sqrt(d2);
	}
	else
	{
		d1 = std::sqrt(d1);
		d2 = 0.5 * q / d1;
	}

	// Set up useful values for the quadratic factors 
	q1 = d1 * d1;
	q2 = -0.25 * c1;
	i = 0;

	// Solve the first quadratic
	p = q1 - 4.0 * (z - d2);
	if (p == 0.0)
		roots[i++] = (float32)(-0.5 * d1 - q2);
	else if (p > 0.0)
	{
		p = std::sqrt(p);
		roots[i++] = (float32)(-0.5 * (d1 + p) + q2);
		roots[i++] = (float32)(-0.5f* (d1 - p) + q2);
	}
	// Solve the second quadratic 
	p = q1 - 4.0 * (z + d2);
	if (p == 0.0)
		roots[i++] = (float32)(0.5f * d1 - q2);
	else if (p > 0.0)
	{
		p = std::sqrt(p);
		roots[i++] = (float32)(0.5f * (d1 + p) + q2);
		roots[i++] = (float32)(0.5f * (d1 - p) + q2);
	}

	if (i > 0)
		x0 = roots[0];
	if (i > 1)
		x1 = roots[1];
	if (i > 2)
		x2 = roots[2];
	if (i > 3)
		x3 = roots[3];

	return i;
}

	} // namespace maths
} // namespace dw

#endif // !DW_MATHS_POLYNOMIAL_INL
